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Abstract 

We present a second-order accurate method to include arbitrary distributions of force densities in the 
lattice Boltzmann formulation of hydrodynamics. Our method may be used to represent singular force 
densities arising either from momentum-conserving internal forces or from external forces which do not 
conserve momentum. We validate our method with several examples involving point forces and find excel- 
lent agreement with analytical results. A minimal model for dilute sedimenting particles is presented using 
the method which promises a substantial gain in computational efficiency. 

PACS numbers: 47. 1 l.-j, 47.15 .G-, 82.70.Dd 
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I. INTRODUCTION 



The numerical integration of the discrete velocity Boltzmann equation provides an efficient 
method for the solution of isothermal, incompressible fluid flows in complex geometries [1]. The 
finite-difference equation generated by the integration scheme is referred to as the lattice Boltz- 
mann equation (LBE). The method can be extended to study multiphase |2j] and multicomponent 
y|] flows, the hydrodynamics of polymers |40 and suspensions |5il6|], and flows under gravity ||7|]. 

In the extensions of the LBE described above, the resulting momentum balance equations con- 
tain additional terms beyond the usual pressure and viscous forces. These represent forces acting 
on the fluid, either from external sources like gravity, or internal sources like the gas-liquid in- 
terface in a two-phase fluid. External sources can add to the total momentum of the fluid, while 
internal sources, being immersed within the fluid, can only exchange momentum with it. Internal 
forces on the fluid, therefore, can always be expressed as the divergence of a stress tensor. This 
formulation encodes the fact that there are no local sources or sinks of momentum. 

The LBE, derived as it is from the Boltzmann equation for a dilute gas, can only faithfully repre- 
sent the hydrodynamics of a fluid with an ideal-gas equation of state and a Newtonian constitutive 
equation ||8|]. One way around this restrictive situation is to use the forced Boltzmann equation 
to represent the additional forces that appear in the extensions described above yj]. So far, this 
idea has been used mainly to model the effects of gravity [70 and gas-liquid interfacial forces in 
non-ideal gases [2]. These correspond to two special types of force distributions: in the case of 
gravity, the force is spatially and temporally constant, while in the case of the gas-liquid interface, 
the force varies both in space and time, but is only evaluated on the nodes of the computational 
grid. The forces in either case are smooth functions of position. 

However, many models of boundaries immersed in fluids require singular distributions of 
forces. Such a description follows, for example, when a gas-liquid interface is described as a 
two-dimensional manifold of zero thickness instead of a three-dimensional volume of space where 
the density changes rapidly. In a similar mathematical idealization, a polymer in a fluid may be 
represented as a one-dimensional curve with a singular distribution of forces |Q] • In yet another 
example, at distances lar ge c ompared to its radius, a sedimenting colloid can be well- approximated 
as a singular point force [j 1 CX ] . Clearly, the range of applications of lattice Boltzmann hydrodynam- 
ics can be greatly expanded if singular force densities, not necessarily located at grid points, can 
be incorporated into the method. 
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In this paper we show how to include force densities having smooth or singular distributions, 
located at arbitrary (in general, off-lattice) points into the lattice-Boltzmann formulation of hy- 
drodynamics. In the following section we first discuss the discrete representation of the forcing 
term in the Boltzmann equation and derive a second-order accurate integration scheme for the 
discrete velocity forced Boltzmann equation using the method of characteristics. In Section UTIl we 
introduce a general distribution of singular forces and using a suitable regularization of the delta 
function obtain a smooth but sharply peaked distribution of forces. The method is exemplified in 
Section Hn] for three common singular force distributions (a Stokeslet, a stresslet, and a rotlet) and 
validated for the Stokeslet case by comparison with fully resolved numerical simulation. Finally, 
we show how our method can be adapted to provide a simplified description of a dilute suspension 
of sedimenting colloids. We end with a summary of our method and discuss potential applications. 



II. MULTIPLE RELAXATION TIME FORCED LATTICE-BOLTZMANN EQUATION 



The LBE may be derived from the Boltzmann equation by a two-step procedure. First, a dis- 
crete velocity Boltzmann equation (DVBE) is obtained by retaining a finite number of terms in the 
Hermite expansion of the Boltzmann equation and evaluating the conserved the moments using a 
Gauss quadrature 111 ill . The discrete velocities {c,} are the nodes of the Gauss-Hermite quadrature. 
This is followed by a discretization in space and time to provide a numerical integration scheme, 
which is commonly called the LBE 012Q. 

Usually a first-order explicit Euler scheme is used to integrate the DVBE, which surprisingly 
enough, gives second-order accurate results [13]. This is so because the discretization error has 
the same structure as the viscous term in the Navier-Stokes equation, whereby it can be absorbed 
by a simple redefinition of the viscosity to give second-order accuracy. The same Euler scheme 
for the forced Boltzmann equation gives a discretization error term which can be absorbed only be 
redefining physical quantities like the momentum and stress [14]. Below, we provide a a straight- 
forward explanation of these redefinitions and show how they are related to the discretization error 
induced by the integration scheme. 

We begin with the discrete velocity Boltzmann equation including an external acceleration field 
F(x,r) 

dtfi + d-Vfi+W-Vcfli = -S?ij(fj-Jf) (1) 
where /, (x, t) is the one-particle distribution function in phase space of coordinates x and velocities 
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C;, J^ij is the collision matrix linearized about the local equilibrium jf and the repeated index j is 
summed over. Mass and momentum conservation require the collision term to satisfy 

t^ij(fj-fj) = 0, E^y(/y-^)c/ = (2) 

i=0 i=0 

while isotropy requires that the S£[j depend only on the angles between c, and Cj yj]. Eq. (0Q) is 
most easily derived by expanding the distribution functions in terms of tensor Hermite polynomi- 
als, truncating the expansion at a certain order, and evaluating the expansion coefficients using a 
Gaussian quadrature II 1511 . lad dimensions, the quadrature is defined by the n discrete velocities c, 
and a set of weights w, giving rise to a DdQn discrete Boltzmann equation ||16|l . Retaining terms 
up to second order in the Hermite expansion is sufficient for isothermal fluid flow problems. The 
equilibrium distribution functions to second-order in the Hermite expansions are 



pv-c,- pvv:Q ; 

P+ c? + 2c? 



(3) 



where the tensor Q ia p = c\ a c^ — c 2 s 8 a p (where Greek indices denote Cartesian directions) and 
c s is the speed of sound. The mass density p and the momentum density pv are moments of the 
distribution function: 



i=0 i=0 



(4) 



To the same order in the Hermite expansion, the discrete representation of the forcing term is given 
by 



[F-V c /]« = -pwi 



F-c; (vF + Fv):Q f 



2c4 



(5) 



Finally, the deviatoric momentum flux tensor 

n 

Saj3= n aj 6 - pc s 8 a p = ^ fiQiap 



(6) 



i=0 



is the second moment of the distribution function. In isothermal models, the higher moments 
represent non-conserved kinetic degrees of freedom, commonly known as ghost modes. In the 
hydrodynamic limit, Eq. (OQ) gives rise to Navier-Stokes behaviour, described by 



p(d,v + v-Vv) = -Vp + T]V 2 v+£V(V-v)+F 



(7) 



where the pressure obeys p = pc^, and the shear viscosity r\ and the bulk viscosity £ are related 
to the eigenvalues of J^ 7 . In practice the algorithm is normally used in a parameter regime where 
the fluid is nearly incompressible (V • v ~ 0). 

To begin our derivation of the numerical scheme we rearrange Eq. (OQ) to obtain 

dtfi + a-Vfi^Rifrt) (8) 

where Ri(x,t) = —J£ij[fj(x,t) — fj(x : t)] +4>;(x,f) represents the effects of both collisions and 
forcing. Eq. © represents a set of first-order hyperbolic equations and can be integrated using the 



method of characteristics fl 1 8h . Integrating over a time interval At we have 

-At 



risr 

fi(x + C{At,t + At) - fi(x, t) = dsRi(x + C{S, t + . 

Jo 



(9) 



The integral above may be approximated to second-order accuracy using the trapezium rule and 
the resulting terms transposed to give a set of implicit equations for the ff. 

fi(x + CiAtj + At) - + c,-Af , t + At) = 

fi{x,t)-^-Ri(x,t)+AtRi(x,t) (10) 

The structure of the above set of equations suggests the introduction of a new set of auxiliary 

nn 

distribution functions yfl, H19I1 

f i (x,t) = f i (x,t)-^-R i (x,t) (11) 



in terms of which the previous set of equations are explicit, 

f,(x + c l At,t + At)=f [ (x,t)+R i (x,t)At. (12) 

This shows that the LBE evolution can be thought of two separate processes: the first is a 
relaxational step in which the distributions fi are relaxed to their "post-collisional" values fi(x, t*), 

fi{x,t*)=f i {x,t)+R i (x,t)At, (13) 

followed by a propagation step in which the post-collisional distributions are propagated along a 
Lagrangian trajectory without further change, 

f 1 (x + c 1 At,t + At)=f l (x,t*). (14) 
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Thus the computational part of the method is most naturally framed in terms of the auxiliary distri- 
butions fi and not the physical distribution functions fi themselves. To obtain the post-collisional 
/; without having to refer to the fi, the latter must be eliminated from Eq. (PT21) . Inverting the 
equations defining the fi in Eq. (fTTT) we obtain 

R, = (1 + ^)7/ [-^ft(A->f) + *y(«,0] • (15) 

Combining this with Eq. (fT2l) we obtain a numerical scheme for the forced discrete Boltzmann 
equation with a general collision operator in terms of the fi: 

i;<x + c i Ar,r + Ar)=/;(x,r) + (l + ^) .H-^(A-A )+4> i (x,0]. (16) 

For a single relaxation time collision operator, where Jz^ ; - = 5y/ t, this takes on a particularly 
simple form 

fi(x + c i At,t + At)=f i (x,t) + -^[-(fi-ff) + T® i (x,t)], (17) 

T + 2 I — 

a result obtained previously by a multiscale expansion of the LBE dynamics Il20n . For a non- 
diagonal collision operator, the collision term is best evaluated in the moment basis. For example, 
using a collision operator in which the ghost modes are projected out and the stress modes relax 
at a rate T _1 , the post-collisional fi {i.e. the RHS of Eq. (fT6l) ) is given by 



fi(x,t*)=wi 



(18) 



where A a , the momentum component of the post-collisional auxiliary distributions, is 



A a = Y,fiC, a +pF a At (19) 

f=0 

and B a p, the stress component, is 

n 

L fiQiap ~ P v aVp + t(v«F j8 + F a Vp ) J (20) 



i=0 



At 



x + At/2 v . =0 

The hydrodynamic variables are moments of the physical distribution fi, but can easily be 
obtained from the auxiliary distributions fi used in the computation, using the transformation rule, 
Eq. (fTTT) . the definitions of the macroscopic variables, Eq. ©, and the constraints of mass and 
momentum conservation, Eq. ©. We obtain 



P = £/i, (21a) 

i=0 

n - At 

pv« = E^c/a + pF a — , (21b) 

;=0 Z 

" At/2 f " - \ 

Sap = J^fiQiap + T + At / 2 I L^'«J3 ~ ^ VaV i6 + T ( v «^ +^ 7 «v j8 ) 1 (21c) 

The equilibria can be reconstructed from p and p v. What appear in the literature as redefinitions 
of momentum and stresses are shown in the above analysis to be discretization errors which vanish 
as A — > tO. This completes the description of the method for the numerical solution of the forced 
LBE. Verberg and Ladd have derived results equivalent to those above using a multiple scale 



analysis of the discrete LBE dynamics [14]. It is not clear to us whether their analysis admits 
singular force densities. However, the above derivation shows that these equations are a reliable 
starting point even in that case. 

The LBE can be extended to situations where the fluctuations in the fluid density and momen- 



tum are important [50 . A consistent discrete kinetic theory of fluctuations was presented in H21H . 
which improves on an earlier algorithm due to Ladd [5], and produces thermodynamically accu- 
rate variances of the local mass and momentum densities. We return to the issue of noise below, 
when we address the representation of Brownian colloids as point particles (Section |TVT). 

IE. SINGULAR FORCE DENSITIES 

In a wide variety of situations, as mentioned in the Introduction, force densities may need to be 
defined off-lattice, and may in addition be singular. Mathematically, such a force density may be 
written as 

F(r) = J f(R)5(r-R)dA (22) 

where the force is localized to some manifold described parametrically as r = R(A) and dX is 
the measure on the manifold. Any numerical method which attempts to deal with such force dis- 
tributions must be reconciled with the singular nature of the force and, for grid-based numerical 
methods, the fact that the position of the manifold need not coincide with the nodes of the grid. 



In a well established numerical method [22], the Dirac delta function in the singular force distri- 



bution is replaced by a regularized delta function which leads to a smooth distribution of forces. 
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Necessarily, this implies that the force is now no longer localized on the manifold but is sharply 
peaked and smooth around it. This smooth force density can now be sampled on the grid using 
the discretized delta function as an interpolant. Thus a representation of Eq. (1221) on the grid is 
obtained from 

F(r)=£f(R fl )5 /, (r-R fl ) (23) 

a 

The crucial ingredient here is the kernel function 8 P , which is a representation of the Dirac delta 
function regularized on the grid. We have followed closely the method described by Peskin I122ll 
where a regularized approximation to the Dirac delta function with compact support is derived: 



f 



(24) 



where h = Ax = Ay = Az is the lattice spacing and f(r) is given by: 



' 3— 2|r| + - N /l+4|r|— 4r 2 



f{r) 



M < l. 



5-2H-V-7 + 12H-^ 1<W< , 



(25) 







\r\ > 2. 



This form is motivated by the need to preserve the fundamental properties of the Dirac delta 
function on the grid (220. A simple closed form approximation to 8 P which is useful for analytical 
work is 



f{r) 



whose Fourier transform is given by: 



i(l + cos(f)) |r|<2, 



\r\ > 2, 



(26) 



f(k) = sinc4^fc + ^sinc (Ank — + ^sinc (Ank + n) . 



(27) 

In this work we combine Eq. (122)) directly with the numerical method described in the previous 
section, giving a well-defined method for incorporating singular and/or off-lattice force densities 
into the lattice Boltzmann hydrodynamics. 



A. Validation 

To validate the method, we compare analytical solutions of the singularly forced Navier-Stokes 
equation against our numerical solutions, using lattice units Ax = Ar = 1 , p = 1 . The most straight- 
forward benchmark is against the initial value problem for the Stokes limit 

d t y = -Vp + T]V 2 v + F(r) (28) 
8 
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FIG. 1: Relaxation of the solenoidal component of the velocity due to a delta function forcing. The simula- 
tions were performed on a 64 3 lattice. Shown are the first 8 (upper curves), the 16th (middle curve) and the 
31st (lower curve) Fourier modes of the velocity field. The solid line is the analytical result, Eq. 



where the nonlinearity has been discarded, incompressibility is assumed, and F(r) = Fq8 (r — Rq) . 
In an infinite system, the solution is obtained in terms of the unsteady Oseen tensor describing the 
diffusion of vorticity Il23ll . In a system with periodic boundary conditions, the Oseen solution 
must be replaced by the Hasimoto solution II24I1 . In contrast to the Oseen solution, the real-space 
Hasimoto solution is not available in a simple closed form but must be evaluated numerically. 
However, the solution in Fourier variables presents no such difficulty, and is in fact identical in 
both cases: 



i 



Tfifc 2 



■(l-kk)-F(k) 



(29) 



Thus, we find it most convenient to compare Fourier modes of the velocity from the numerical 
solution against the solution above. In particular, this provides a neat way to evaluate the perfor- 
mance of the method at different length scales. In FIG. CD we compare the numerical data (points) 
to the theoretical result for a regularized force monopole using the approximation to the Peskin 
delta function, Eq. dTTT) (solid line). 

The results show excellent agreement with the theoretical curve for low k modes, where we 
expect the momentum to behave hydrodynamically. The departure from hydrodynamic behaviour 
increases progressively with the wavenumber, as expected from previous studies on the hydro- 
dynamic behaviour of the LBE [25Q. However, there is a significant range of length scales over 
which our model reproduces hydrodynamic behaviour, which is not less than the scale over which 
hydrodynamic behaviour is obtained in the unforced LBE Ii26ll . 
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(b) 

FIG. 2: (Colour online) Velocity around a symmetric point-force dipole, normalized by a characteristic 

Fa 



speed for that distance from the dipole, v c = ^f-p. ■ FIG. |2(a)] shows velocity along lines parallel to the 
forces, at several separations. Points are simulation results, lines theoretical predictions with no free param- 
eters. In figure |2(b)[ the upper half shows the simulated velocity field. The lower half shows isosurfaces 
of the magnitude of the velocity difference between simulation and theory at values of 25% and 50%. The 
colouring (online) depends upon the magnitude of the difference field and is shown as a percentage of v c in 
the colour bar. The force dipole is oriented vertically and positioned in the centre of the volume. 

By combining elementary monopoles, discrete representations of higher multipoles can be gen- 
erated. For example, the discrete Stokes doublet, a dipole of two point forces, can be constructed 
out of monopoles of magnitude F and separation a and is often used as a simplified representa- 
tion of a neutrally buoyant, steadily moving self-propelled particle. In figures |2(a)| and \2(b)\ we 
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FIG. 3: (Colour online) Velocity field around a regularized rotlet, normalized by v c (see FIG. El)- Upper 
half: simulated velocity field. Lower half: isosurfaces of the magnitude of the velocity difference between 
simulation and theory. Isosurfaces are at values of 12.5%, 25% and 37.5%. The colouring (online) depends 
upon the magnitude of the difference field and is shown as a percentage of v c in the colour bar. The rotlet is 
oriented with the forces in a horizontal plane and positioned in the centre of the volume. 

compare the velocity response of such a dipole to theoretical predictions, finding good agreement 
away from the immediate vicinity of the forces. In FIG. [3l we show a velocity field plot for the 
antisymmetric force dipole, or rotlet, which may be used as a representation of an object which 
rotates due to an external torque. This requires the use of four, rather than two, point forces; we ar- 
range these in a swastika-like fashion (whose axes can be aligned in an arbitrary direction without 
significantly affecting the flow produced). This cancels a spurious stresslet component that arises 
from our regularization of the 8 function for any dipole in which the forces are not collinear with 
the separation vector. 

The above examples show that the regularized delta function provides an useful way of incor- 
porating arbitrary distributions of singular forces into the lattice Boltzmann method, capable of 
dealing with internal as well as external forcing. 

IV. A STOKESLET MODEL FOR DILUTE COLLOIDS 



The dynamics in a dilute sedimenting suspension, despite a century of investigation, still 
presents open questions 111 411 . The problem, even for a hard-sphere suspension, is unusually dif- 
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ficult due to the long-ranged, many-body nature of the hydrodynamic interaction. Moreover, the 
flow can develop structural features at large length scales, and the role of inertia, while usually 



negligible at the particle scale, may be significant at those scales 112711 . The Stokes approximation 
of globally vanishing Reynolds number cannot thus be justified a priori in a sedimenting sus- 
pension. The full hydrodynamic problem including inertia for both fluid and particles was first 
simulated by Ladd using a novel lattice Boltzmann method [28]. This method, though possibly 
the most competitive for fully resolved particles, remains computationally expensive. A consid- 
erable simplification of the hydrodynamics is possible if only the lowest order multipole of the 
force distribution induced on the particle surface by the no-slip boundary condition is retained. 
This principle was exploited previously to develop representations of polymers as strin gs of po int 



particles which were then coupled to an LB fluid [4 , 



2911 . A similar idea has been used [30, 



aid to 



represent resolved colloids with a mesh of point particles covering their surfaces. However in the 
current work we simplify further, treating each colloid as a single point particle (thereby sacrific- 
ing all near-field effects) .In the colloidal context, this model was first introduced by Saffman 
the finite sized particles are replaced by a singular force monopole, the Stokeslet, located at the 
nominal centre of the particle. In Saffman's original model, both the fluid and the particles have 
no inertia. In keeping with the comments above, our model retains inertia for the fluid, while ne- 
glecting it for the particle and for hydrodynamics at the particle scale. We thus have a momentum 
balance equation, 

p^v = -Vp + T]V 2 v + £F 1 .5(r-R 5 ) (30) 

s 

where the sum includes contributions from the s = 1 . . .N particles located R ? and acted upon 
by external forces F s . In the absence of particle inertia accelerations vanish, and the particle 



coordinates are updated directly using the first Faxen relation 113211 



R s = v~(R,) + -^-, (31) 

which relates the centre of mass velocity of the particle R^ to the external force on it F 5 . The 
background velocity v co (R iS ) is the fluid velocity at the location R ? in the absence of the 5-th 
particle. The above two equations provide a complete specification of a model of sedimenting 
spheres, valid in the dilute limit, for dynamics at long wavelengths. 

The lattice Boltzmann implementation of this model proceeds by first replacing the Dirac 
delta function with the regularized delta functions to obtain a force density at the grid points 
F(r) = £ ? F. v <5 p (r — R. v ). Since the LBE evolves the total fluid velocity v(r) due to all particles, 
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the background fluid velocity v°° must be obtained by a careful subtraction procedure. In the ab- 
sence of fluid inertia at the particle scale this can be accomplished as follows. By definition, the 
fluid velocity at a node v(r) is the sum of the background velocity at the node v°°(r) and the ve- 
locity due to the 5-th Stokeslet located at R$, v(r) = v°°(r) +v iS (r,R lS ). The background velocity 
field at the location of the particle can be obtained using the same interpolation kernel as used for 
the force, v°°(R J ) = J^ r v°°(r)5 /> (R s — r), and using the previous relation can be written as 

v°°(R s ) =v(R s ) -£v s (r,R s )S p (R s -r) (32) 

r 

Appealing only to linearity and dimensional analysis, the sum above can be expressed as 

W'^-^SSS) <33) 

In Appendix [A] we derive this result and show that the lattice parameter depends only on the 
system size L and on the form of the regularization and interpolation kernels; it is independent of 
viscosity r\ and of the radius a. Using Eq. (|33~T) . the update equation for the Stokeslet positions 
can now be expressed in terms of the interpolated fluid velocity, without any reference to the 
background velocity, 

R, = v(R,) + -4rrl ( 34 > 



6nri \a az,(Ry) 

Notice that replacing the background velocity in the Faxen relation with the actual fluid velocity 
induces an effective backflow, leading to a renormalized hydrodynamic radius, 

1 = 1-1 (35) 

The numerics thus places a constraint a C flL on the allowed values of the hydrodynamic radius 
a. This numerical constraint encodes the condition that the grid points must be in the far-field of 
the Stokeslet, the limit in which the background velocity can be obtained from the fluid velocity 
by subtracting a monopole contribution. In our simulations, we operate well within this limit. 

This almost completes the description of the lattice Boltzmann implementation of our Stokeslet 
model of sedimenting particles. The only free parameter is the hydrodynamic radius a of the 
particles, which decides how fast they sediment for a given force F s . As shown below, the lattice 
parameter cll can be calculated analytically as a function of system size. We find it convenient 
to fit it using a procedure described in Appendix lAl Finally, to address Brownian motion of our 



colloids, we need to use the FLBE of bill which imparts an appropriate thermal noise spectrum 



to the fluid. Because of the renormalization of a, the resulting diffusivity is generally not correct 
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FIG. 4: (Colour online) Response of a single particle to an impulsive force P in a periodic box. Main figure: 
response shown for a range of viscosities (points) compared to the theoretical prediction v = n p(nvt)- 3 l 2 at 



long times 1230 . The inset shows the effect of varying the size of the simulation box. Deviations from the 
prediction become significant at approximately 250, 1000 and 4000 timesteps for box sizes of 32, 64 and 
128, respectively. This is consistent with the expected scaling, X ~ L 2 /rj. 

unless a further noise term is added that is the counterpart of the correction. The details are 
explained in Appendix [B] 



A. Benchmarks 



Our first benchmark addresses the dynamics of a single impulsively started particle, without 
noise. From unsteady hydrodynamics, we know that the asymptotic decay of the particle velocity 
varies as t~ d / 2 in d dimensions [23]. In FIG. |4] we display the decay of the particle velocity, for 
a single hydrodynamic radius (0.05 LU), but several values of the fluid viscosity. In all cases, we 
see the correct asymptotic behaviour, until the particle is begins to interact with its image, due to 
the periodic boundary conditions. (This interpretation is supported by the scaling of the time x at 
which the deviations become significant: we find x ~ jr as shown in the inset to FIG. 0]) In other 
words, our particle model correctly captures the low and intermediate frequency behaviour of the 
particle mobility, but cannot capture the high-frequency behaviour correctly, since that depends 
on the way vorticity diffuses in the immediate neighbourhood of the particle, a regime which is 
excluded in our model. 

Our next benchmark involves collective motion of a set of particles, and thus directly probes 
the hydrodynamic interaction between particles. In FIG. [5] the mean sedimentation velocity of a 
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FIG. 5 : (Colour online) Steady sedimentation velocity of a simple cubic array of spheres, normalized by the 
Stokes sedimentation velocity of a single particle, vq. The separation, L, is expressed in terms of the fitted 



particle radius, a. The solid line is the theoretical result v = vq/(1 + &£), with k = 2.84 1124 
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(a)r = 0? St (b)/ = 1650f St (c)f = 3300f St 

FIG. 6: Crowley instability of a 2-dimensional lattice of sedimenting particles. Images are generated in the 
comoving frame and particle size has been exaggerated for illustrative puiposes. 



periodic array of spheres is shown, as a function of volume fraction. There is excellent agreement 
with the theoretical result of 113 311 . The model is also able faithfully to capture instabilities due to 
collective hydrodynamic flow. In FIG. [6] the instability of a falling 2D lattice of spheres in three 
dimensions B4I1 is captured, at least qualitatively, by our model. 

For most problems, all Reynolds numbers below some (situation-dependent but small) value 
give rise to equivalent behaviour, as discussed in detail in previous work 113 511 . Following protocols 
discussed there, we have compared the normalized velocity field u = \/vq (with vo the sedimen- 
tation velocity of an isolated colloid) for a number of simulations of a single sedimenting sphere 
with periodic boundary conditions (see figure d) in order to explore the range of Reynolds num- 
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(a) Re = 1(T 6 (b) Re = 1(T 4 (c) Re = 1(T 2 

FIG. 7: Contour plots of the normalized velocity field u for a reference simulation at very low Re = 10~ 6 
(top); and velocity difference fields for Re = 10~ 4 (middle) and 10~ 2 (bottom). These are for a point- 
like colloid sedimenting in a 32 3 box with periodic boundary conditions. Reference case: contour interval 
0.02vo. Re = 10~ 4 case: contour interval 5 x 10~ 6 . Re = 10~ 2 case: contour interval 5 x 10~ 3 . 

ber at which our algorithm gives acceptably accurate results. Our 'reference' simulation has a 
very small Re = 10~ 6 such that we can be confident it is the in the Stokesian limit 113 511 . This 



is shown in the panel |7(a)| Panels |7(b)| and |7(c)| show the normalized velocity difference fields 



between the reference case and simulations with Re = 10~ 4 and Re = 10 -2 , respectively. In the 
simulation with Re = 10~ 4 , the magnitude of the difference is everywhere less than 2 x 10~ 5 , a 
negligibly small error. In the simulation with Re = 10~ 2 , we find |Au| < 0.01 throughout the bulk 
of the domain; only in the immediate vicinity of the particle does it become larger. This sug- 
gests that this Reynolds number is sufficiently low to give 'realistic', although not 'fully realistic', 



behaviour H35H . Since reaching very low Reynolds number requires paying a larger cost in compu- 
tational time, and there are other sources of percent-level error in the code, Re = 10 2 is probably 
a reasonable compromise between accuracy and run-time, for studies in the low Reynolds number 
limit. 



B. Comparision to a fully resolved LB algorithm 



As a final benchmark, we have compared the behaviour of our sedimenting particle model with 
a fully resolved colloid simulation code using the algorithm of Nguyen and Ladd 113 611 . (For full 
implementation details see 0711 .) At dilute concentrations, the paths of the particles are almost 
indistinguishable between the two simulations when plotted graphically. This is shown in figure [8] 
for volume fraction = 3.0 x 10~ 3 ; note that the largest differences occur when the density of 
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FIG. 8: (Colour online) Both visualizations show the particles at their position in the point-particle simu- 
lation after 10 Stokes times. Left: lines show the trajectories from the starting configuration. Right: lines 
show the difference between fully-resolved and point like algorithm. Parameters for the two systems are 



shown in TablelU See also Movie 1 [38]. 



particles is large locally, when the implicit assumption of our model that the particles are always at 
separations large compared to their radius is no longer valid. We cannot expect both simulations to 
give the same trajectories for long times, since the small differences between algorithms will cause 
exponential separation of trajectories owing to the positive Lyapunov exponent of the system. 
However, from a plot of the mean difference in position between the two simulations against time, 
we can see excellent agreement for several Stokes times, and until at least ten Stokes times for 
sufficiently dilute systems (FIG. [9]). 



V. DISCUSSION 



The focus of this work has been to derive and validate a general method for addressing singular 
forcing in the LBE, with specific application to the simulation of point-like particles. We have 
shown the method to agree well with analytic results, where available, and with fully-resolved 
particle algorithms at low concentrations and Reynolds number. Additionally, due to its careful 
construction, the regularized 5-function provides a good interpolation scheme, minimizing veloc- 
ity fluctuations as the particle moves relative to the computational grid. Indeed we find that for 
sedimenting colloids the trajectories are much smoother in our Stokeslet algorithm than for the 
fully resolved simulation. In the latter, the discretization renders particles hydrodynamically as- 
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time / ts 



FIG. 9: (Colour online) Main figure: mean absolute difference in position between fully-resolved and 
point-like sedimenting particles. Inset: mean absolute difference in velocity between the two simulations. 
Parameters for the two systems are shown in Table HI 



Parameter 


Fully resolved 


Point particle 


Lattice size L 


96 


9 


Particle radius a 


1.25 


0.117 


Viscosity 77 


1 

6 


1 

6 


Density p 


1.0 


1.0 


Reynolds number Re 


0.01 


0.01 


Stokes time is 


937 


8 


Number of particles 


321 


321 


Volume fraction 


3.0 x 10~ 3 


3.0 x 10~ 3 



TABLE I: Parameters for simulations used to compare fully-resolved and point-like algorithms; see text and 
figures [8] and [9] 



pherical with shapes that vary as they move across the lattice Ii37ll . Absence of such irregularities 
in the Stokeslet code may make this generally preferable at small volume fractions. 

For dilute suspensions of sedimenting colloids, our new algorithm can thus perform simulations 
of accuracy comparable to (or even better than) that of a fully resolved code, but at vastly reduced 
computational cost. As shown in Table U similar particle numbers and volume fractions can be 
simulated with an LB lattice that is smaller in linear dimension by a factor A ~ 10. (This is the 
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ratio of the particle radii in the two simulations.) The computational time to update the particle 
positions is essentially negligible, so that the CPU time needed to perform one LB time step is 

2 

decreased by a factor of A 3 ; moreover the Stokes time, T5 = scales as A 2 . The latter sets the 
time basic time scale for evolution of sedimentation trajectories, so that for this problem we expect 
a speed-up of 0{X 5 ) ~ 10 5 . This should allow us to study the sedimentation behaviour of dilute 
systems with tens of millions of particles; we hope to pursue this avenue elsewhere. 
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APPENDIX A: SUBTRACTION PROCEDURE 



We derive here the correction factor that arises from replacing the background velocity with the 
fluid velocity in the Faxen relation. The velocity at a node due a regularized Stokeslet located at 

Rs 1S k r 

v,(r,R,) = 1 £ e -ls p (k;R s )(l -kk) -F, (Al) 
This velocity interpolated from the neighbouring nodes to the location of the particle is 

£v>,R,)5>-R,)=F,. £ ^^^(l-kk)^(r-R,) (A2) 
Completing the spatial sum, we get the for the interpolated Stokeslet velocity, 

X>.,(r, RX(r - «.) = F, • £ (1 - tt) = (A3, 

which shows that the offset parameter ai obeys 



a L (R s ) L 3 k ^ 



is indeed independent of viscosity, and particle radius a, but depends on the lattice size L and the 
numerical implementation of the regularization and interpolation. 
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APPENDIX B: NOISE 



Considering a spherical sedimenting particle in the Langevin picture, we can write 

mR= -6^T]a(R-v°°) +F+C(0 (Bl) 
and taking the inertialess limit gives 

ft = ^ R > + 6^<" <B2) 

We note that in this equation, noise only comes in through the Gaussian random variable £ and 
that in particular the fluid velocity v TC is completely deterministic. 
The update rule for our model particle is 

R = v(R) - + (B3) 

6nria L 6nr\a 

which is sufficient for the infinite Peclet number regime. If one uses fluctuating LB, then the 
interpolated velocity v(R) contains a noise component. However, we do not expect the magnitude 
of the noise to be appropriate for a particle of the desired radius a, since the random component 
of the velocity has no dependence on the radius of the particle. In fact, the variance of this noise 
is that expected for a Brownian particle with the same radius as the offset parameter (> a) and 
we use this fact to determine its value from a diffusion "experiment" on an unforced particle, as 
explained below. 

Knowing that the particle will otherwise diffuse as one with a much larger radius, we add a 
white noise term to the update rule for the model 

R = v(R) - — ?— + + f'(0 (B4) 
6nr\a L 6Ttr\a 

The variance of the extra noise is determined, by the requirement of satisfying the fluctuation- 
dissipation theorem, to be 

<#(0Cj(O> = (---) Sij8(t-t'). (B5) 

To determine the value of the offset parameter a^, we set up a simulation of a single unforced 
particle in periodic boundary conditions at finite temperature and disable the extra noise term 
discussed above, giving 

R = v(R) (B6) 
20 



as the equation of motion of the particle. 

We let the simulation equilibrate for the characteristic time for momentum to diffuse across the 
box size, T L 2 /ri, before recording the displacement as a function of time. This is repeated for 
a number of different starting positions relative to the LB grid and a plot of (r 2 ) vs. t is used to 
estimate the diffusivity. We then use the Stokes-Einstein relation to derive a radius and use this as 
the offset parameter. We then test to ensure that this gives the correct sedimentation behaviour of 
a particle. 
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